function fig3_13

%  uses box method to solve
%  y'' + p(x)y' + q(x)y= f(x)   for xL < x < xR

% set boundary conditions
xL=0; yL=1;
xR=1; yR=2;

nx=10000;
ep=0.0001

a0=0; b0=0;

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.05,'P',0.04)

x=linspace(0,1,nx);
h=x(2)-x(1);
% calculate the coefficients of finite difference equation
a=zeros(1,nx-2); b=zeros(1,nx-2); c=zeros(1,nx-2);
for i=1:nx-2
    xm=x(i+1)-0.5*h;
    xM=xm+h;
    a(i)= -2 - 0.5*h*(p(xM,a0,b0,ep)-p(xm,a0,b0,ep)) + 0.25*h*h*(q(xM,a0,b0,ep)+q(xm,a0,b0,ep));
    b(i)= 1 - 0.5*h*p(xm,a0,b0,ep) + 0.25*h*h*q(xm,a0,b0,ep);
    c(i)= 1 + 0.5*h*p(xM,a0,b0,ep) + 0.25*h*h*q(xM,a0,b0,ep);
    f(i)= 0.5*h*h*(rhs(xm,a0,b0,ep) + rhs(xM,a0,b0,ep));
end
f(1)=f(1)-yL*b(1);
f(nx-2)=f(nx-2)-yR*c(nx-2);

% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);
y=[yL, y, yR];

% interior layer solution
%x3=1/3-20*sqrt(ep);  x33=1/3+20*sqrt(ep);
x3=0.2;  x33=0.6;
nxx=200;
xx=linspace(x3,x33,nxx);
al=1;  ar=2^(10/9)*exp(1/3);
alp0=0.5*(ar+al)*gamma(4/9)/ ( sqrt(pi)* (6*ep*exp(2))^(1/18) );
bet0=sqrt(1.5/pi) * (ar - al)*gamma(17/18)/ ( (6*ep*exp(2))^(1/18) );
for i=1:nxx
    xxx=(xx(i)-1/3)/sqrt(ep);
    yy(i) = alp0*M(1/18,1/2,-1.5*xxx^2)+bet0*xxx*M(5/9,3/2,-1.5*xxx^2);
end

plot(x,y,'-','Linewidth',1.2)
hold on
plot(xx,yy,'--r','Linewidth',1.4)
box on
grid on
axis([0 1 0.9 4])

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('Solution','Interpreter','latex','FontSize',14)

legend(' Numerical Solution\,\,',' Interior Approx','Location','NorthEast','FontSize',13,'Interpreter','latex')

%exportgraphics(gcf,'/Users/mark/Desktop/figk2.pdf')


function g=q(x,a,b,ep)
g=x/ep;

function g=p(x,a,b,ep)
g=(3*x-1)/ep;

function g=rhs(x,a,b,ep)
g=0;


% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
    y(j) = y(j) - v(j)*y(j+1);
end




function g=M(a0,b0,x0)

% this code will not work if b0 is a negative integer
%if (b0==round(b0)) & (b0<0)
%	error('The code doesn't work if b is a negative integer');
%	return;
%end

% use identity:  M(a,b,x)=exp(x)*M(b-a,b,-x)
if x0< 0
    x=-x0; a=b0-a0; b=b0;
else
    x=x0; a=a0; b=b0;
end

xmax=20;

% evaluate M(a,b,x) for x ge 0

%  use exact expressions for M if possible
if (x==0) | (a==0)
    g=1;
elseif a==b
    g=exp(x);
elseif (a==1) & (b==2)
    g=(exp(x)-1)/x;

    %  using the complete poly isn't the smartest idea if a is a very large
    %  negative integer but its done anyway
elseif (a==round(a)) & (a<0)
    N=-a;
    g=1;
    cs=1;
    for in=1:N
        cs=cs*(a+in-1)/((b+in-1)*in);
        g=g+cs*x^in;
    end

    %  otherwise, for small x use series definition
elseif x<xmax
    N=8000;
    g=1;
    cs=1;
    for in=1:N
        cs=cs*(a+in-1)*x/((b+in-1)*in);
        g=g+cs;
        if abs(cs/g)<1.0e-12
            break
        end
    end

    %  and use asym approx for large x
else
    N=400;
    g0=1;
    cs=1;
    for in=1:N
        cs=cs*(b-a+in-1)*(1-a+in-1)/(in*x);
        g0=g0+cs;
        if abs(cs/g0)<1.0e-12
            break
        end
    end
    q1=gamma(b)/gamma(a);
    g=x^(a-b)*exp(x)*q1*g0;
end

%  adjust answer if x0<0

if (x0<0)
    g=g*exp(x0);
end





